VAN_DER_POL

Overview

The VAN_DER_POL function numerically solves the Van der Pol oscillator, a classic non-conservative dynamical system with nonlinear damping. Named after Dutch physicist Balthasar van der Pol, who discovered it while studying vacuum tube circuits at Philips in the 1920s, this oscillator is fundamental to the study of nonlinear dynamics and chaos theory.

The Van der Pol oscillator is governed by the second-order differential equation:

\frac{d^2x}{dt^2} - \mu(1 - x^2)\frac{dx}{dt} + x = 0

where x is the position coordinate, t is time, and \mu is the nonlinearity and damping parameter. This equation exhibits self-sustaining oscillations: near the origin the system is unstable (energy is added), while far from the origin damping removes energy. This balance produces a stable limit cycle to which all initial conditions converge when \mu > 0.

The function reformulates this as a two-dimensional system suitable for numerical integration:

\begin{aligned}
\frac{dx}{dt} &= y \\
\frac{dy}{dt} &= \mu(1 - x^2)y - x
\end{aligned}

The behavior depends strongly on \mu:

When \mu = 0, the system reduces to a simple harmonic oscillator
For small \mu, oscillations are nearly sinusoidal
For large \mu, the system exhibits relaxation oscillations with slow buildup and rapid release phases

This implementation uses SciPy’s solve_ivp function from the SciPy library (GitHub) for numerical integration. Multiple solver methods are available including RK45 (default explicit Runge-Kutta), Radau, and BDF for stiff problems. For more theoretical background, see the Van der Pol oscillator article on Wikipedia or the detailed treatment on Scholarpedia.

The Van der Pol oscillator has applications beyond electrical engineering, including modeling biological rhythms such as cardiac pacemakers and neural action potentials (through the related FitzHugh-Nagumo model), seismology, and phonation studies.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=VAN_DER_POL(x_initial, y_initial, mu, t_start, t_end, timesteps, solve_ivp_method)

x_initial (float, required): Initial position x of the oscillator.
y_initial (float, required): Initial velocity y of the oscillator.
mu (float, required): Nonlinearity and damping parameter (mu >= 0).
t_start (float, required): Start time of integration.
t_end (float, required): End time of integration.
timesteps (int, optional, default: 10): Number of timesteps to evaluate.
solve_ivp_method (str, optional, default: “RK45”): Integration method for solve_ivp.

Returns (list[list]): 2D list with header row [‘t’, ‘X’, ‘Y’]. Each subsequent row contains time and variable values. str: Error message if input is invalid or integration fails.

Examples

Example 1: Demo case 1

Inputs:

x_initial	y_initial	mu	t_start	t_end	timesteps
2	0	1	0	10	10

Excel formula:

=VAN_DER_POL(2, 0, 1, 0, 10, 10)

Expected output:

t	X	Y
0	2	0
1.111	1.418	-0.8404
2.222	-0.136	-2.304
3.333	-2.03	-0.0125
4.444	-1.459	0.8189
5.556	0.03088	2.206
6.667	2.007	0.06472
7.778	1.452	-0.8189
8.889	-0.04559	-2.218
10	-2.007	-0.0507

Example 2: Demo case 2

Inputs:

x_initial	y_initial	mu	t_start	t_end	timesteps	solve_ivp_method
2	0	1	0	10	10	RK23

Excel formula:

=VAN_DER_POL(2, 0, 1, 0, 10, 10, "RK23")

Expected output:

t	X	Y
0	2	0
1.111	1.418	-0.8404
2.222	-0.134	-2.302
3.333	-2.007	0.02457
4.444	-1.417	0.8428
5.556	0.1394	2.307
6.667	2.007	-0.02921
7.778	1.415	-0.8442
8.889	-0.145	-2.313
10	-2.007	0.0339

Example 3: Demo case 3

Inputs:

x_initial	y_initial	mu	t_start	t_end	timesteps
2	0	2	0	10	10

Excel formula:

=VAN_DER_POL(2, 0, 2, 0, 10, 10)

Expected output:

t	X	Y
0	2	0
1.111	1.652	-0.4272
2.222	0.9881	-0.9084
3.333	-1.58	-2.796
4.444	-1.839	0.3576
5.556	-1.337	0.5947
6.667	-0.02533	2.569
7.778	1.988	-0.2417
8.889	1.581	-0.4576
10	0.8292	-1.106

Example 4: Demo case 4

Inputs:

x_initial	y_initial	mu	t_start	t_end	timesteps
1	1	0.5	0	1	10

Excel formula:

=VAN_DER_POL(1, 1, 0.5, 0, 1, 10)

Expected output:

t	X	Y
0	1	1
0.1111	1.104	0.8773
0.2222	1.194	0.7353
0.3333	1.267	0.58
0.4444	1.323	0.4169
0.5556	1.36	0.2513
0.6667	1.379	0.08803
0.7778	1.38	-0.06835
0.8889	1.364	-0.2144
1	1.333	-0.3515

Python Code

from scipy.integrate import solve_ivp as scipy_solve_ivp
import numpy as np

def van_der_pol(x_initial, y_initial, mu, t_start, t_end, timesteps=10, solve_ivp_method='RK45'):
    """
    Numerically solves the Van der Pol oscillator system of ordinary differential equations.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        x_initial (float): Initial position x of the oscillator.
        y_initial (float): Initial velocity y of the oscillator.
        mu (float): Nonlinearity and damping parameter (mu >= 0).
        t_start (float): Start time of integration.
        t_end (float): End time of integration.
        timesteps (int, optional): Number of timesteps to evaluate. Default is 10.
        solve_ivp_method (str, optional): Integration method for solve_ivp. Valid options: RK45, RK23, DOP853, Radau, BDF, LSODA. Default is 'RK45'.

    Returns:
        list[list]: 2D list with header row ['t', 'X', 'Y']. Each subsequent row contains time and variable values. str: Error message if input is invalid or integration fails.
    """
    # Validate input types
    try:
        x0 = float(x_initial)
        y0 = float(y_initial)
        mu_val = float(mu)
        t0 = float(t_start)
        t1 = float(t_end)
    except (TypeError, ValueError):
        return "Invalid input: All initial values, parameters, and times must be numeric."

    try:
        ntp = int(timesteps)
    except (TypeError, ValueError):
        return "Invalid input: timesteps must be an integer."

    # Validate input values
    if t1 <= t0:
        return "Invalid input: t_end must be greater than t_start."
    if ntp <= 0:
        return "Invalid input: timesteps must be a positive integer."

    valid_methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']
    if solve_ivp_method not in valid_methods:
        return f"Invalid input: solve_ivp_method must be one of {', '.join(valid_methods)}."
    # Create time array for evaluation
    t_eval = np.linspace(t0, t1, ntp)

    # Van der Pol ODE system
    def vdp_ode(t, state):
        """Van der Pol oscillator equations.

        Args:
            t: Time (not used but required by solve_ivp).
            state: [x, y] state vector.

        Returns:
            [dx/dt, dy/dt] derivative vector.
        """
        x, y_var = state
        dxdt = y_var
        dydt = mu_val * (1 - x**2) * y_var - x
        return [dxdt, dydt]
    # Integrate
    try:
        sol = scipy_solve_ivp(
            vdp_ode,
            [t0, t1],
            [x0, y0],
            method=solve_ivp_method,
            dense_output=False,
            t_eval=t_eval
        )
    except Exception as e:
        return f"Integration error: {str(e)}"

    if not sol.success:
        return f"Integration failed: {sol.message}"
    # Format output: header row then data rows
    result = [["t", "X", "Y"]]
    for i in range(len(sol.t)):
        t_val = float(sol.t[i])
        x_val = float(sol.y[0][i])
        y_val = float(sol.y[1][i])

        # Check for NaN or inf in output
        if not np.isfinite(t_val) or not np.isfinite(x_val) or not np.isfinite(y_val):
            return "Invalid output: NaN or inf detected in solution."

        result.append([t_val, x_val, y_val])

    return result

Overview

Excel Usage

Examples

Python Code

Online Calculator